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ABSTRACT 

We report on simulations of the formation of the first stars in the Universe, 
where we identify regions of hot atomic gas (/h 2 < 1CT 6 ) at densities above 
10~ 14 g cm -3 , heated to temperatures ranging between 3000 and 8000 K. Within 
this temperature range atomic hydrogen is unable to cool effectively. We describe 
the kinetic and thermal characteristics of these regions and investigate their ori- 
gin. We find that these regions, while small in total mass fraction of the cloud, 
may be dynamically important over the accretion timescale for the central clump 
in the cloud, particularly as a chemical, rather than radiative, mechanism for 
clearing the polar regions of the accretion disk of material and terminating ac- 
cretion along these directions. These inherently three-dimensional effects stress 
the need for multi-dimensional calculations of protostellar accretion for reliable 
predictions of the masses of the very first stars. 

Subject headings: galaxies: formation; stars: formation; ISM: H II regions; cos- 
mology: theory 



1. Introduction 



As the physical model that governs simulations of metal-free stars expands to include 
more chemical and radiative processes, the collapsing regions in which they are expected to 



form have revealed a variety of interesting phenomena (jYoshida et al.ll2006t iTurk et al.ll2008 



Clark et al. 


2010; 


Turk et al. 


2009; 


Stacv et al. 


2010) 



metal- free halos are strongly governed by the re actions that govern the formation and disso- 



ciation of molecular hydrogen at high densities ( ITurk et al.ll2010bl submitted). At densities 
below 10 -16 g cm -3 , the formation of mol e cular hydrogen is gov erned by electron-catalyzed 



association of H and H( lAbel et al.l 119971 ; Idover fc Abel! 120081 ). However, at densities of 
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10 16 g cm 3 and above, the formation of molecular hydrogen is dominated by three-body 
reactions, primarily where the third body is atomic hydrogen: 

H + 2H H + H 2 
H + H 2 -> H + 2H. 

The reaction rates for these t wo reactions have been shown not only to be uncertain to an 
order of magnitude or more ( iGloverl 120081 ) but these uncertainties may have a substantial 



impac t on the structure and mass-distributions of collapsing primordial clouds (ITurk et al. 



2010b|). Every molecule formed through this reaction releases 4.48 eV in thermal energy into 



the gas, and because the rate of formation decreases with increasing temperature, the overall 



several decades in density ( 


RiDamonti & Abel 


2004 




Omukai & Nishi 


1998 




Galli & Palla 


1998 




Palla et al. 


1983 




Glover 


2008 


). In ID, spherically symmetric calculations it is in- 



evitable that all high-density gas undergoes this transition. However, in full 3D calcula- 
tions, the mechanism of transformation is more complex, owing to variations in the density- 
temperature distribution and the velocity structure of the collapsing cloud. 

We present results of high-resolution computer simulations of metal-free star forming 
regions. In these simulations we identify regions of high-density (> 10~ 15 g cm~ 3 ) atomic gas 
whose presence may contribute to a chemical, rather than radiative, mechanism for slowin g 
or terminating accretion along the poles of primordial accretion disks ( IMcKee fc Tanll2008l ). 



2. Simulations 



We report on a single simulation drawn from a suite of calculations, designed to sample a 
variety of formation environments. We centered the simulation volume on the location of the 
peak density of the first dark matter halo of mass 10 6 M Q to form and generated three nested 
subg rid levels, fo r an in i tial effective resolut i on of 1024 3 . Our initial conditions are generated 

e exception that we hay e 



as m 



Abel et all (I2002f ); lO'Shea fc Norman! ( 120071 ). with the notab 



updated our cosmo logical parameters to WMAP 7 (IJarosik et al. 



Thes e simulations were c onducted using the code Enzo (lO'Shea et al. 



2010] lTuAet~al1l2010ah. 



2004 



Bryan fc Norman 



19981 ; iBryan et al.l 120011 ) , usi ng a 9-specie s primordial chemistry model imp l emen ting the 
TW0STEP method described in lVerwerl ( 119941 ) and in the appendix to lTurk et al.l ( 120091 ). While 
we focus on the effects of three-body molecular hydrogen formation with atomic hydrogen 
acting as the third body, we also solve the set of rate equations where molecular hydrogen 
acts as the third body. At densities of approximately 10~ 15 — 10~ 13 g cm -3 , primordial gas 
becomes optically thick to the cooling of molecular hydrogen fro m ro- vibrational transitions . 
In our simulations, this is applied using the prescription given in lRipamonti &: Abel! ( 120041 ). 
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(Turk et al. 


2010b: 


Sobolev 


1960 



Yoshida et al.ll2006l ). The cooling radiation is reduced by: 

£iines,thick(T) = Li ineSj thin x rnin(l, (n/n ) 
with n = 8 x 10 9 amu cm" 3 and (3 = -0.45 ( Ripamonti fc AbeJl200~4l ). 



We discuss a relatively slow-collapsing halo that we can examine with some depth 



but this halo is n ot extraordinary, nor does it fragment (ITurk et al.ll2009t IStacy et al.l 12010 



Clark et al.ll2010l ). Furthermore, the phenomena discussed in this paper are present in a suit e 
of simulations, and their frequency will be discussed in the forthcoming ITurk et al.l (j2010al ). 
The final output of the calculation at z — 17.4, has peak density of 2.0 x 10~ 9 g cm" 3 at 
the 29th level of refinement. The entire simulation has 3.3 x 10 7 computational elements, 
including 2.3 x 10 6 with cell size less than 35 AU. 



We also introduce an improved refinement criterion as compared to previous simulations. 
We require the numerical resolution to be one sixteenth of the local Jeans length. However, 
the calculation of the Jeans length is done assuming the gas is at T m j n = 200 K even when 
it has not cooled so low. As a consequence we have (T/200K) 9//2 more resolution elements 
in the material before it forms the cold phase. We refer to this as the cold Jeans refinement 
criterion, which for all the material above the minimum temperature is more stringent then 
the standard criterion. This modification ensures high resolution at high densities even 
before the material forms the cooling cores. In addition, we also refine on overdensity in 
both baryonic and dark matter content at a factor of 4.0 which ensures high resolution in 
the lower density material that forms first star hosting mini halos. 



3. Results 

As the gas collapses and is converted to molecular hydrogen, we see a spread in the 
temperature and the molecular hydrogen fraction as a function of density, as shown in 
the upper left and upper right panels of Figure [Q At densities between approximately 
10~ 15 g cm -3 and 10 -10 g cm -3 the mass is distributed over a broad range of temperatures, 
while the corresponding spread in molecular hydrogen fractions is an order of magnitude or 
more. This results in a broad variation of the rate coefficients governing the transition to 
molecular hydrogen. For higher temperature gas, the overall transition occurs more slowly, 
resulting in less efficient ro-vibrational cooling. In some regions of phase space, the onset of 
the transformation of the gas is delayed until the onset of optical thickness. 

In the upper left panel of Figure [1] we plot the distribution of mass as a function of 
density and temperature. The clear trend is for the gas to follow standard tracks through 
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Fig. 1. — We show the distributions of the thermodynamic and kinetic states of the inner 
10 4 AU of a collapsing primordial Population III star-forming region gas as a function of 
different variables: mass distribution as a function of density and temperature (upper left); 
mass distribution as a function of the declination (co-latitude) angle from the local angu- 
lar momentum vector and temperature (upper middle); mass- weighted average molecular 
hydrogen mass fraction as a function of density and temperature (upper right); mass dis- 
tribution as a function of radius from the center and an entropy-like quantity (lower left); 
mass distribution as a function of the radius and the local radial velocity of the gas (lower 
middle). In the lower right panel, we show a slice through the center of the simulation of the 
molecular hydrogen mass fraction. The field of view in the slice image is 1000 AU. Colorbars 
showing the molecular hydrogen fraction in each of the molecular hydrogen fraction plots 
and the mass in each pixel of the mass-distribution plots have been placed on the far right, 
in the upper and lower panels respectively. 
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Fig. 2. — Slices at three different epochs (columns) in the collapse of the primordial cloud. 
The field of view in all images is 2000 AU on a side. The axes of alignment of the three rows 
have been chosen in an orthogonal frame (not the computational axes) where the normal 
vector of the top image is aligned with the angular momentum vector of the innermost 
2000 AU of the disk, and the middle and lower rows are mutually orthogonal. Overplotted 
are velocity vectors and a contour at T = 1000 K. From left, the times of the images are 
3415, 1700, 1000, 320, and 85 years prior to the final output of the simulation. 
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phase s pace, similar to those sho wn in lAbel et al.l ( 120021 ) ; lYoshida et al.l ( 120061 ) ; iTurk et al. 
( 120091 ); 10 'Shea &: Norman! (120071 ) . However, we note the strong presence of low-mass regions 



of high-temperature gas, particularly at densities between 10 15 g cm 3 and 10 



13 — 3 

g cm . 



In the upper middle panel of Figure [T] we have plotted the distribution of mass in the 
innermost 10 4 AU of the halo as a function of the temperature and the cosine of declina- 
tion from the average angular momentum vector at that radius. We have chosen cosine to 
ensure equal latitudinal area across the entire domain. In the normalization used here, the 
two extremes ( — 1.0 and 1.0) correspond to the polar regions and the middle region (0.0) 
corresponds to the disk region. All gas above a temperature of roughly 2000 K is confined 
within 60 degrees from the poles defined by the local angular momentum vector. 

The upper right panel of Figure [T] shows the mass- weighted average molecular hydrogen 
mass fraction as a function of density and temperature. All of the gas above 2000 K is mostly 
or completely dissociated and free of molecules. The lower-temperature gas follows a typical 
transition from atomic to molecular states. 

In the lower left panel of Figure [TJ we have plotted the distribution of mass with respect 
to radius from the densest point of our cloud (x-axis) and an entropy-like quantity (S, y-axis) 
defined by 



s _ h T 

where T is the temperature, \x is the mean molecular weight and p is the density. 
Between 10 AU and 100 AU the high-temperature regions are visible as a broadening in the 
entropy profile, indicating a difference in the entropy gradient within the collapsing region. 

In the lower middle panel of Figure [1] we plot the mass distribution as a function of 
distance from the densest zone and the radial velocity with respect to that zone. The 
velocity of the central zone in the calculation has been subtracted prior to calculation of the 
radial velocity. Several independent flow lines are visible, including several with net outward 
motion, particularly within ~ 50 AU. The bulk flow is inward and the spread in velocities 
very close to the central point is expected as a result of both rotation and incoherent flow. 

The lower right panel of Figure [1] gives a slice through the core, aligned with the y-axis, 
1000 AU on a side, where the color corresponds to the molecular hydrogen mass fraction. 
We note the obvious molecular hydrogen-dominated disk and atomic high temperature polar 
regions. 

In Figure El we display slices through the primordial cloud that are aligned with a new 
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orthogonal reference frame such that the normal to the image plane is aligned with the 
angular momentum vector of the innermost 2000 AU of the region. Overlaid are velocity 
vectors and contours at temperatures of 1000 K (dotted) and 2000 K (solid). Each column 
corresponds to a different epoch in the cloud's collapse, where time of 0.0 years is the final 
output of the calculation. The upper panel shows disk settling at densities of approximately 
10~ 14 g cm -3 . This disk undergoes little evolution over the displayed epochs in the top, 
face-on view, In the middle and lower panels, gas is flowing along the poles of the disk 
to the central point. High-temperature gas develops in these polar regions, at and above 
2000 K, although we note that the largest pockets of hot gas are at 1000 K and extend nearly 
1000 AU. 



4. Discussion 

4.1. Chemical Origin and Stability 

The formation of molecular hydrogen at densities greater than 10 -16 g cm -3 is domi- 
nated by the three-body formation mechanism, where two hydrogen atoms and a catalyst 
produce a single molecule of hydrogen. Formation through this channel produces an excited 
molecule of hydrogen, which rapidly collisionally de-excites to its ground state, depositing 
4.48 eV of thermal energy in the gas. Conversely, collisional dissociation of a molecule of 
hydrogen removes from the gas 4.48 eV of thermal energy. As the formation and destruction 
rates of molecular hydrogen via these processes decrease and increase, respectively, with the 
temperature of the gas, this reduction in thermal energy carries with it a corresponding re- 
duction in the rate at which molecular hydrogen associates. This process can be viewed as an 
increase in the specific heat of the gas; every molecule of hydrogen that is to be dissociated 
requires the introduction of 4.48 eV of energy. The equilibrium density of atomic hydrogen, 
at a fixed temperature and subject only to the three-body process and its inverse, can be 
written as: 

™H,eq = 4^(^13 - ^/8fci 3 fc22TlH,H2 + ■ i 1 ) 

where A43 is the dissociation rate in units of cm 3 s _1 and k 22 is the association rate in 
cm 3 s -1 , nH, eq is the atomic hydrogen density in cm -3 and nn,H 2 is the numbe r dens ity of the 
mixe d hydrogen gas. In our simulation, we used values for k 22 from iGloverl (120081 ) and fci 3 



from iMartin et al.l (119961 ). As shown in the upper right panel of Figured! when gas exceeds 
roughly 2300 K at a density of 10~ 14 g cm -3 , the destruction rate dominates and the gas 
rapidly dissociates, but reaching or breaching that temperature barrier requires an influx of 
energy to dissociate any molecules of hydrogen. We note that while the temperature at which 
dissociation dominates will be weakly dependent on density, as seen in Equation [TJ Having 
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less energy "locked up" in hydrogen molecules lowers the energy barrier to dissociation, as 
the gas can shed less energy via dissociative cooling. 

In the absence of shocks, asymmetries or other perturbations of the gas, as in spherically- 
symmetric, ID calculations, molecular hydrogen that is allowed to cool will simply slow its 
collapse if its temperature reaches the temperature at which the dissociative process domi- 
nates. The gas will loiter until radiative processes are able to shed enough thermal energy to 
allow it to collapse further. However, in our simulations, the g as may be s u bject t o asymmet 



(2009 


); 


Clark et al. 


(2010 


); 


Stacy et al. 


(2010 





Abe 



et al. 



(J2002); iTurk et al. 



physical model includes optical thickness to ro-vibrational cooling and our simulation re- 
solves inhomogeneities of the density and temperature of the gas, both factors that change 
this evolution. 

With a smaller specific heat as a result of a lower molecular fraction, the molecular 
hydrogen can be dissociated by minor shocks and other perturbations from infalling gas; 
once it has reached high temperatures, the energy barrier to associating molecular hydrogen 
is just as high, but with no mechanism for shedding thermal energy until compressional 
heating drives the temperature to ~ 10 4 K. In the time-series plots in Figure [2] we show 
that the flow, while largely quiescent, features strong asymmetric inflow as well as fast 
infall along the poles. This produces a shocked, high-temperature gas which is subsequently 
compressionally heated. 

At this point, the gas will be able to cool efficiently enough that it will likely not 
heat further, but it will also be unable to form substantial amounts of molecular hydrogen. 
This transition is shown in the upper right panel of Figure HJ where the high-temperature, 
shocked regions in the upper left panel are shown to have fully dissociated their reservoir 
of molecular hydrogen, and thus the only available coolant. We also draw attention to the 
power-law stratification in the molecular hydrogen fraction as a function of density and 
temperature, indicating an equilibrium molecular hydrogen fraction at those temperatures. 

In the bottom left panel of Figure [Tj we can identify two features in the entropy of the 
gas. The first is that the entropy overall decreases with decreasing radii; the cloud is globally 
stable against convection, which would require a reversal of the entropy gradient. However, 
the warm gas behaves adiabatically, retaining essentially constant entropy with decreasing 
radius. During subsequent accretion onto the protostellar cloud, this gas will continue to be 
confined by colder, radiatively efficient gas. Although globally stable against convection, at 
later times these adiabatic regions may undergo convective instability. In the final output 
plotted in Figure |2] (bottom rightmost image) we note that a pocket of gas of ~ 1000 K has 
become distended and nearly disconnected from its innermost component; this may be the 
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start of a buoyant path, up and to the right in the image plane. 



4.2. Radial Velocity and Disk Structure 

In the lower middle panel of Figure [TJ we can identify several distinct "tracks" along 
which gas is flowing. Each of the hot, dissociated regions is characterized by amplification 
of the inward radial velocity and a relatively small radial extent. Gas is falling inward onto 
the hot atomic bubbles, where it dissociates, increases in pressure, and then builds up on 
the inward side of the velocity spikes. This variance in the radial velocity of the shocked, 
dissociated regions results in gas that is collapsing less slowly; in the comoving frame of its 
surroundings, it may even be rising, although we cannot determine this at the stopping point 
in our simulation. 

As seen in Figure HJ in this simulation a clear disk is formed, accentuated by the strong 
differences in the chemical and thermal states of the gas in the disk and in the polar regions. 
Examining the angular distribution of gas (Figure [lj upper middle panel) we see that the 
shocks are essentially confined to within or near the polar regions of the disk, demonstrated as 
well in the lower right panel of Figure [TJ The mechanism by which Population III protostellar 
feedback will break out of the disk, as well as the density structure of a settled disk both 
require an understanding of this material, as the gas flowing along the poles will necessarily 
interact with these high-entropy regions. In past theoretical works discussing the accretion 
mechanism onto primordial protostars, radiative breakout has been identified as the sole 
mechanism by which accreting material flowing along the poles could be slowed or stopped 



(IMcKee fc Tanll2008l ). This scenario may change if the polar regions are largely composed of 
high-temperature, atomic gas that is buoyant or unstable to convection. We propose that a 
chemical component may contribute to changes in the accretion flow even before the onset 
of radiative feedback from the central protostar, at times later than those probed in this 
simulation. The decreased molecular opacities in the polar region will allow protostellar 
accretion luminosity to reach larger radii, feeding back on the accretion flow. 



5. Conclusions 

The hydrodynamics of disk formation around a protostar lead to high temperature 
atomic polar regions. Until compressional heating drives the temperature of these hot, 
dissociated regions to 10 4 K, they have no opportunity for radiative cooling. The initial 
formation of these bubbles is through shocks dominated by large scale flows, and the chemical 
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processes that govern their evolution will continue to enable their formation. We see these 
dissociated regions prior to the formation of the central protostar, and we extrapolate that 
they will continue to persist until the protostar begins to accrete. Clearing material from the 
polar regions of the accretion disk, previously thought to be exclusively through radiative 
feedback mechanisms, is necessary to reverse accretion flows along the polar directions and 
enable the breakout of ionizing radiation. However, the process of terminating accretion 
along the polar regions of the accretion disk may be assisted by these hot regions, either by 
the dredging of material via buoyant or convective rising of material, or simply through the 
chemical and kinetic processes that prevent the region from associating molecular hydrogen. 
This may contribute to a much earlier breako ut of radiation from the protostellar core than 
found in the models of iMcKee fc Tan! (120081 ). although (as noted in that paper) the total 
mass of gas flowing along the poles is much lower than that accreted from the disk. 

The spread in the chemical, thermal and kinetic states of the gas as it regimes such 
as where we see these shocking regions, must be resolved to understand the mechanism by 
which accretion disks develop and matter is processed. One-dimensional codes, or codes that 
assume a universal equation of state as a function of density, are unable to do so. Further- 
more, simulation codes that are unable to resolve the development and persistence of shocks, 
or that are otherwise unable to apply adequate resolution, are equally unable to resolve this 
formation and settling of a disk and its subsequent e volution. Future simulations that bypas s 
the Courant condition, such as sink particles (as in Idark et al.l feoioh : Istacv et all (bold )) 
while still resolving the relevant chemical and hydrodynamical processes, will be necessary 
to study the long term impact of these hot, atomic polar regions. 
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